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Abstract — For a power system operating in the vicinity of the 
power transfer limit of its transmission system, effect of stochastic 
fluctuations of power loads can become critical as a sufficiently 
strong such fluctuation may activate voltage instability and lead 
to a large scale collapse of the system. Considering the effect of 
these stochastic fluctuations near a codimension 1 saddle-node 
bifurcation, we explicitly calculate the autocorrelation function 
of the state vector and show how its behavior explains the 
phenomenon of critical slowing-down often observed for power 
systems on the threshold of blackout. We also estimate the 
collapse probability/mean clearing time for the power system and 
construct a new indicator function signaling the proximity to a 
large scale collapse. The new indicator function is easy to estimate 
in real time using PMU data feeds as well as SCADA information 
about fluctuations of power load on the nodes of the power grid. 
We discuss control strategies leading to the minimization of the 
collapse probability. 

Index Terms — Blackout prevention, emergency control, phasor 
measurements, power system stability, voltage stability, wide-area 
measurements and control. 



I. Introduction 

IT is well known that small stochastic fluctuations of power 
load, although usually negligible in the vicinity of a stable 
operating point, may potentially lead to a large scale cascading 
failure if the power system operates close to a saddle-node bi- 
furcation point JTJ, (2). Early detection and mitigation of such 
failures is a problem of utmost importance in contemporary 
world of constantly increasing power demand, where power 
grids often operate in a precritical regime. 

Statistical properties of aggregated power load, their influ- 
ence on static and dynamic characteristics of power systems 
remain the subject of extensive studies for the last three 
decades, see [1J for the review. A very considerable atten- 
tion has been given to stability analysis of power networks 
based on Lyapunov theory of dynamical systems, as it was 
recognized early that the proximity of an operating point 
to saddle-node and/or Hopf bifurcations of power systems 
signals about approach to an instability of the base state (3). 
In particular, saddle-node bifurcations were associated to the 
phenomenon of voltage collapse as well as loss of synchronism 
H). Correspondingly, a multitude of stability criteria based on 
the estimation of global Lyapunov functions of power systems 
have been introduced (see for example classical papers 0, 
[6 1, |7]), many corresponding indicator functions being already 
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used in industry for preventive control and dynamic security 
assessment of power grids (T]. Despite its popularity, energy 
function analysis of power systems is not entirely free of 
drawbacks: (a) typically, a knowledge of the values of system 
variables on all nodes of the network is needed in order to 
estimate the global Lyapunov function in a given operating 
regime, and (b) even if such tremendous amount of data 
is available in real time, the estimation process has a high 
computational cost, which becomes especially critical in the 
proximity of a large scale collapse. 

Assuming in the present work that the operating point of the 
power system under consideration is close to a saddle-node 
bifurcation, we extend the standard approach to Lyapunov 
stability analysis of power systems, taking into account that 
for a typical power grid without any specific structural sym- 
metries the center manifold is one-dimensional, i.e., saddle- 
node bifurcation has codimension 1. This observation allows 
us to explicitly calculate the autocorrelation function of system 
variables in the operating regime near the bifurcation point 
and find an approximate expression for the mean clearing 
time/probability of a large scale failure of the power system. 
Using the latter, we construct a new indicator function of 
proximity to the voltage collapse/loss of synchronism, which 
is significantly easier to estimate in real time than the global 
Lyapunov function, especially if the system operator receives 
real time PMU data as well as SCADA information about 
fluctuations of power demand on individual nodes. 

This manuscript is organized as follows. In Sec.[ll]we review 
the structure-preserving model used in the present paper to 
describe dynamic behavior of system variables in the vicinity 
of the saddle-node bifurcation point. The autocorrelation func- 
tion of system variables is explicitly calculated in the Sec. [HI] 
while a new estimate for the mean clearing time is given in 
the Sec. IV In the Sec. [V] we perform the validation of our 
model by checking the derived formulae against numerical 
simulations of the IEEE 39 bus (New England) power system 
described in J8). Finally, the Sec. [Vljis devoted to discussion 
of control strategies for minimizing the collapse probability 
and conclusions. 



II. The model 

To describe dynamics of system variables in the vicinity of 
a saddle node bifurcation point, we use the system of coupled 
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swing equations on (P, V) nodes of the grid 
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and power flow equations on (P, Q) nodes 
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where as usual 8i is a voltage phase on a bus i, Hi is an inertia 
constant for a generator on the node i, a.i denote frequency 
controls on the (P, V) nodes with generators, Vi is a voltage 
magnitude on a bus i (for the (P, V) nodes, Vi = Ei). Real and 
reactive power loads on (P, Q) nodes are generally (weakly 
changing) functions of frequency 9i, voltage Vi and voltage 
change rate V. If only dynamics of system variables at the 
vicinity of a stable operating point is to be considered, it is 
sufficient to keep first leading orders in the Taylor expansions 
of real and reactive power loads in powers of their arguments. 
This corresponds to the structure preserving model of ||9l 
extended to the case of load dependence on the voltage change 
rate Vi. 

The loads have a fixed power factor fc^. They fluctuate with 
time, and fluctuations of loads on different nodes of the grid 
are statistically independent. We assume that their correlation 
properties are Gaussian: 



(6P i (t)SP j (t')) = B i (t-t , )6 i 
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One simple example of the function P,; is exponential P; (t — 
t') = Ai exp(— \t — t'\/ri), where the characteristic time scale 
t is the decay time of temporal correlations of power load 
fluctuations. 

Since the load power factor is fixed, fluctuations of reactive 
loads are related to ([4]) as 

(SQiWSQjffl) = k ^ L B i (t-t')6 ij , 

so that only fluctuations of 5 Pi are needed to be considered. 

Note that Q describes the autocorrelation function of a sta- 
tionary process. The property of stationarity (or translational 
invariance in time) t — > t+5t, t' — > t'+St is only approximate, 
as behavior of loads and their random fluctuations features 
natural cycles (for example, day/night). However, we are only 
interested to study a short time scale (tens of seconds and 
minutes) behavior of the state vector, when Q is perfectly 
applicable. 

After solving the power flow equations, finding the base 
state and linearizing equations ([TJ, |2]), Q about the stable 
operating point, one finds the matrix stochastic differential 
equations (SDE) 



Mx + Vx + JCx = Ax = SP. 



(5) 



Here x is the system state vector including voltage phases 
and magnitudes, the matrix A4 describes inertial properties 
of generators connected to the grid, the matrix T> — primary 



frequency controls on (P, V) nodes as well as frequency and 
V dependence of power loads on (P, Q) nodes, and finally JC 
is the power flow Jacobian encoding all static properties of 
the power system. The system |5]) of SDE will be the main 
subject of our study. 

III. Calculation of autocorrelation function 

Generally, the autocorrelation function of the system vector 
x can be found as 

du> 
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where 



(X(U)X T {-LU)) = ^- 1 (W)S(CJ)(^( W ))- 



and the Fourier-transformed system matrix is 

A(u) = -Muj 2 + jVuj + JC. 

The correlation properties of fluctuating loads are given by the 
matrix 

B{uS) = (SP(ui)SP i '(w)) = B(ui)l. 

The value of the integral |6]l is determined by the singularities 
of the integrand in the complex ui plane, which in particular 
include zeros of detA(ui) and det A(—uj), as well as the 
singularities of B{ui). In the regime of interest, dynamics 
of system variables is largely dictated by very slow modes 
(see below). Power loads fluctuate rapidly compared to the 
time scales of this slow dynamics, their correlation in time 
described by Q becomes negligible, and one can simply write 
B(oS) = B ■ 1, where P^ = L°° dtBij(t). For the case 
Bij(t) = Ai exp(—t/Ti)5ij considered in the previous Section, 
one expect that the time scales Tj are the shortest in the system, 
and one effectively has P^- = A,-/ri<Jy. 

Naturally, the singularity closest to the real oj axis de- 
termines behavior of the autocorrelation function |6| at late 
times. To identify it, we note that at a saddle-node bifurcation 
point several eigenvalues of the power flow Jacobian JC vanish. 
There exists only one such vanishing eigenvalue e — > if the 
center manifold of the power system is one-dimensional [12|. 
Performing the eigenvalue decomposition of the inverse power 
flow Jacobian 

ic-i = J2kKW, 

i 

where Aj are eigenvalues of JC and a; and bi are corresponding 
left and right eigenvectors, we see that if the operating point 
of the power system is close to the saddle-node bifurcation 
point, JC^ 1 is approximately given by 



JC' 1 = -ba 1 
e 



(7) 



where e is the eigenvalue of JC vanishing at the bifurcation 
point. Dots denote the contribution of all the other eigenvalues 
of JC, which is at most of the order 0(1) in powers of e 
(although for large systems it can be numerically large). We 
are particularly interested in the situation when the 2-norm 
1 1 -6a — /C _1 ||2 <C ||/C _1 ||2 and the asymptotic representation 
holds well. 
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In the case under consideration, the leading singularity of 
the integrand in (j6j) coincides with a zero of detA(u>) (or 
det A(— u>) depending on the sign of the time difference t—t'). 
Such singularity is a simple pole by assumption that the center 
manifold of the power system is one-dimensional. Solving the 
equation det*A(w) = det(— A4uj 2 + iDuj + JC) = by means 
of perturbation theory in the small parameter e, one finds that 
the mode determining behavior of the autocorrelation function 
at late times is given by 
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The frequency of the leading mode is purely imaginary^] 

Estimating the integral (|6]l in the vicinity of this leading 
singularity, one finally finds (to the leading order in e) 



(x(0)x T (t)) 
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Note that in this operating regime near the codimension 1 
saddle-node bifurcation point dynamics of system variables 
x(t) is directly related to static characteristics of the power 
system such as the lowest eigenvalue e of the power flow 
Jacobian. 

As follows from |9]), the right eigenvector b of the power 
flow Jacobian JC determines the preferred direction in the 
phase space of a power system near a codimension 1 saddle- 
node bifurcation fl2l . ifTTIl . In turn, the left eigenvector a 
shows fluctuations of the loads on which nodes mostly de- 
termine stochastic dynamics of the state vector: the nodes i 
with larger components m are the ones, where fluctuations of 
the power loads SPi, SQi influence the dynamics of the state 
vector stronger. 

As the operating point of the power system approaches the 
power transfer limit and the lowest eigenvalue e of the power 
flow Jacobian gets smaller, the rate of decay of the autocorrela- 
tion function |9} decreases quickly. This observation is known 
as the phenomenon of critical slowing-down observed during 
large scale failures in power grids [13], lfT4l . 

Finally, we would also like to emphasize that fluctuations 
of the leading mode, the component of the system vector x 
along b direction, become strongly amplified as e decreases, 
since the amplitude of the leading mode near the bifurcation is 
proportional to e -1 / 2 . One of the main goals of system control 
in this regime is to decrease the amplitude of the leading mode 
by adjusting the values of the matrix element a T T>b (and, if 
possible, a T Ba). This can be done using primary frequency 
control on (P, V) nodes as well as utilizing resources of load 
following on (P,Q) nodes as we discuss below. 

IV. Mean clearing time and probability of 

VOLTAGE COLLAPSE 

Taking the results of the previous Section into account, 
it is also possible to explicitly calculate the mean clearing 
time for an operating point clos^] to the bifurcation point. 

' An oscillating contribution disappears from the mode for sufficiently small 
e < (a T Vb) 2 /4a T Mb. 

2 But not too close as will become clear from the subsequent discussion. 



Namely, considering a scalar reduced system variable z(t) 
defined according to 

z(t) = b T x(t), 



one finds the following equation of motion for it: 



a T Mb ■ z + a T Vb ■ z + ez + a T Ybb - 2 
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where the matrix JTyfe = -g^- is the first derivative of 
the power flow Jacobian w.r.t. the system variables and . . . 



denote higher derivatives of fC. The equation (Hi is a scalar 
SDE describing the process of activation of an over-damped 
unharmonic oscillator. Therefore, one can simply apply the 
classical result by Kramers |[T8l to calculate the probability of 
collapse/the mean clearing time: 
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Z{a T YW?(a T Ba) , 

Estimates similar to this one were previously found using 
Lyuapunov stability analysis of power systems (see for exam- 
ple [6), 0). Yet, the expression (12i is different from these 



well-known results in one important respect: the mean clearing 



time (12 1 depends only on a small number of parameters 
of the power system, unlike the global Lyuapunov function 
which is a functional of dynamical characteristics of the power 
system. This is so because the number of degrees of freedom 
of a power system gets enormously reduced in the vicinity of 
bifurcation, and all the information about the power system in 
the vicinity of the bifurcation point is aggregated through a 
small number of quantities. These quantities (such as matrix 
elements a T Vb and a T Ba or the lowest eigenvalue e of the 



power flow Jacobian) entering (12i can be straightforwardly 



estimated from the static state analysis using the measurements 
of the correlation function (x(0)x T (<)) of the state vector 
(containing the information about the vector b, the matrix 
element a T T>b and e) and local SCADA data of power 
consumption on individual nodes of the grid (which encode 
information about the matrix B) and therefore do not require 
any a priori knowledge of the power system structure. Such 
measurements will make it easier to estimate the mean clearing 
time for the power system in real time using the expression 

We also emphasize that the result ( [12} is strictly speaking 
applicable only to the case of power systems with one- 
dimensional center manifold, and an estimate similar to ( p"2j ) 
would be impossible to make if the dimension of center man- 
ifold of the power system is 2 or larger. What is necessary for 
the derivation of ( [T2| is the equilibration law known in physics 
as the principle of detailed balance, and it does not apply 
for system matrices of generic power grids. However, because 
of tremendous reduction of the total number of degrees of 
freedom in the vicinity of a codimension 1 bifurcation point, 
the power system becomes effectively one-dimensional (there 
is only one relevant degree of freedom), and the principle of 
detailed balance is always valid for one-dimensional systems. 

V. Numerical simulation of IEEE 39 power system 

We shall now perform the verification of the model de- 
scribed in the Section [II] comparing the expression (|9]l for the 
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correlation function with results of numerical simulations of 
the IEEE 39 (New England) power system. The strategy is as 
follows. First, the saddle-node bifurcation point is localized 
using continuation power flow procedure [15| implemented 
in PSAT Toolbox for Matlab |16|, then it is identified more 
precisely using MATPOWER 4.1 library for Matlab [171- 
An operating point close to the bifurcation is chosen, static 
power flow equations — solved using MATPOWER, the power 
flow Jacobian — calculated, and finally stochastic differential 
equations of the Section [TT] are solved using the implicit Euler 
algorithm. The number N r of realizations of the random load 
fluctuations is 1000 (we consider only a stationary regime 
where averaging the system variables x(t) over realizations 
of the statistical ensemble Q of random loads is expected to 
be equivalent to averaging of x(t) over time). 

The load parameter A = 2.12 is chosen (if A = 1 
corresponds to the operating point described in |8|). This 
corresponds to a 0.9% displacement in terms of the generated 
real power as compared to the critical regime^] As usual, 
the bus 31 is a slack bus. The state vector includes voltage 
phases on the nodes 1-30, 32-39 and voltage magnitudes — 
on the nodes 1-29. The same A4 matrix chosen as presented 
in (8), with the base frequency /o = 60 Hz. The T) matrix is 
diagonal with components d = 1 p.u. corresponding to voltage 
phases on (P, Q) buses, d = 0.1 p.u. corresponding to voltage 
magnitudes on (P, Q) buses and d = 10 p.u. — to phases on 
(P, V) buses. The loads are allowed to fluctuate only on the 
buses 3, 10 and 21, with the same characteristic amplitudes of 
fluctuations \J~B = 0.1 p.u. For our choice of A the smallest 
eigenvalue of the power flow Jacobian JC is e « 0.1811 p.u, 
with the ratio of2-norms | lAC^ 1 — e^^al | 2 /| lAC^ 1 1 1 2 « 0.0717. 
Therefore, the expansion in powers of e is feasible, albeit e 
is not too small by itself (only smallness of dimensionless 
quantities dictates where the perturbation theory is applicable). 

The numerical results for the autocorrelation function 
(z(t)z(tf)) of the reduced system variable z(t) = b T x{t) 
as well as the scalar autocorrelation function (x T (i)x(t')) 
are presented on the Fig. [T] The full correlation function 
of the system vector (x T {t)x(t')) clearly follows the one of 
reduced variable (z(t)z(t')) except the initial short interval 
of time, where contributions from sub-leading modes into 
the autocorrelation function are not small. This confirms our 
expectation that fluctuations of system variables orthogonal 
to the direction of the vector b become suppressed in the 
vicinity of the saddle-node bifurcation. Note that for any given 
finite number of realizations N r the autocorrelation function 
(x T (i)x(i')) fluctuates stronger than (z(t)z(t')) because it 
also accounts for fluctuations of the system vector orthogonal 
to the vector b, although in the limit of infinite number 
of realizations both correlation functions coincide with each 
other. 

For our choice of parameters the value of the amplitude 
of the correlation function (x T (t)x(t')) is a T Ba/2ea Tr Db rj 
5.96 • 10~ 4 , while the best fit value recovered from [I] is 
5.26 • 10~ 4 . The autocorrelation decay time r = a T T>b/e ss 



17.3 sec also coincides well with the best fit value r w 17.0 
sec and an integral estimate for the average correlation time 

/*' dt{x T {t)x{t'))/{x T {t)x{t)) « 18.1 sec. 



- (x'(t)x(f)) 

- (Z(t)2(f)> 




Et/a T Db 



Figure 1. Relation between the logarithms of the autocorrelation function 
{z(t)z(t')) of the reduced system variable b T x(t) (red) and the full autocor- 
relation function (x T (t)x(t')) (blue). The initial moment of time t is chosen 
to be t = a T T>b/e. 



VI. Discussion and path forward 

In the present paper we discuss a general power grid 
without any particular structural symmetries operating near the 
power transfer limit of its transmission system. A saddle-node 
bifurcation present in the phase space of such power system 
typically has codimension 1, which implies that dynamics 
of state variables near such bifurcation point is essentially 
determined by a single degree of freedom which we denote as 
reduced system variable. We consider how stochastic fluctua- 
tions of power loads exactly influence this dynamics, explicitly 
calculating (a) the autocorrelation function of the state vector 
Q and (b) the probability of a large scale failure of the power 
system in the vicinity of the saddle-node bifurcation point/the 
mean clearing time for the system 



12 1. Although such an 



operating regime was extensively studied in literature for the 
last 30 years (see (T) and references therein), both results (a) 
and (b) are new to our knowledge. Also, using (a) we are able 
to quantify the phenomenon of critical slowing-down recently 
discovered in power grids operating on the threshold of a large 
scale cascading failure fl3l . fl4l . 

Using the expression for the mean clearing time ([12]), one 
can introduce a new simple indicator function ( fT3] > signaling 
proximity of a power system to a large scale failure. Indeed, 
as follows from the expression (12 1, collapse probability is 



» 1 



(13) 



The total generated real power at A cr i t 
13.67 GW. 



2.13569843 is approximately 



exponentially small when the expression 

e 3 a T Vb 
c ~ (a T Tbb) 2 a T Ba 

and vise versa. Even if the operating point of the power 
system under consideration is close to the power transfer limit 
(i.e., e nearly vanishes), the mean clearing time could still 
be exponentially large when the matrix element a T Vb is large 
and/or the matrix element a T Ba — small. This in turn implies 
that the probability of a large scale failure is low. 
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As follows from ( |13) , a strategy for a system operator to 
keep the power system close to the power transfer limit of 
its transmission system while minimizing the probability of a 
large scale collapse would be to maximally utilize resources of 
load following [21 1 (effectively reducing the matrix element 
a T Ba) and/or adjust, if it is possible, parameters of the primary 
frequency control on (P, V) nodes (increasing the value of the 
matrix element a T T>b). Note that a T T>b cannot be increased 
indefinitely: such increase for example would imply a growth 
of the correlation time r = a T T>b/e leading in turn to possible 
issues related to interference between primary and secondary 
frequency controls. 

There exist several reasons why the operating regime of a 
power system discussed in the paper seems rather attractive 
despite the fact that the operating point is located near the 
threshold of system instability. First of all, in such a regime, 
the power grid and its transmission system are utilized with the 
best effectiveness possible: the throughput of the transmission 
system is maximal, nearly all energy which the grid is capable 
to produce is consumed (assuming the minimal power losses 
in the transmission system). Therefore, such operating regime 
of the grid is optimal from the economical point of view given 
the system remains under control and the operating costs of 
control systems are not too high. 

Second, synchronism of the power system actually holds 
rather well in the regime under consideration. Indeed, the 
degree of synchronism is determined by the correlation func- 
tion of the operating frequency (8uji(t)8ijjj{t')), which for the 
nodes i with large representation in the vector b is given to 
the leading order in e by 



ebi(a T Ba)bj 
(a T Vb) 3 



exp 



e\t-t'\ 
a T Vb 



(14) 



Control countermeasures necessary to maximize the value 
of the indicator function ( 13 1 and minimize the collapse 



probability are the ones which also maximize the value of the 
matrix element a T T>b (and/or minimize a T Ba) and therefore 



decrease the amplitude of the correlation function (14i. 

Finally, since in the regime under consideration dynamics of 
the power system is essentially described by a single degree 
of freedom, the complicated problem of state recovery can 
be effectively solved by using real time data feeds from 
synchrophasor measurement units as well as a much larger 
(and much slower acquired) volume of SCADA data about 
fluctuating aggregated power loads. The same applies to the 
problem of system identification [22]: as we have discussed, 
measurements of the autocorrelation function (x(0)x T (t)) 
using PMU data feeds directly provide information about the 
matrix elements a T T>b and a T Ba, the lowest eigenvalue of the 
power flow Jacobian, corresponding left and right eigenvectors 
a and b. Therefore, such measurements can potentially allow 
for the approximate recovery of the system matrix in the 
vicinity of the bifurcation point or at least its part responsible 
for the dynamics of the leading mode of the state vector. We 
shall discuss this observation more extensively elsewhere. 
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